#include "DoubleMatrix.h"

static const double PI = 3.1415926535897932384626433832795;
static const double aa = 0.5;
static const double bb = sqrt(0.5)*cos(PI/8);
static const double cc = sqrt(0.5)*cos(3*PI/8);

static const double DCT4x4[16]=
{
		aa,		aa,		aa,		aa,
		bb,		cc,		-cc,	-bb,
		aa,		-aa,	-aa,	aa,
		cc,		-bb,	bb,		-cc,
};

DoubleMatrix DoubleMatrix::A(4, 4, DCT4x4);
DoubleMatrix DoubleMatrix::At = A.Transposed();

DoubleMatrix::DoubleMatrix()
{
}

DoubleMatrix::DoubleMatrix(unsigned int newHeight, unsigned int newWidth, const double* data): Matrix(newHeight, newWidth, data)
{
}

DoubleMatrix::DoubleMatrix(const Matrix<double>& rhs)
{
	Matrix::operator=(rhs);
}

DoubleMatrix& DoubleMatrix::ForwardDCT()
{
	if (GetWidth() != 4 || GetHeight() != 4)
		throw "DoubleMatrix: Invalid DCT dimension";

	// Y = A * X * At;
	DoubleMatrix result = A;
	result *= *this;
	result *= At;

	*this = result;

	return *this;
}

DoubleMatrix& DoubleMatrix::InverseDCT()
{
	if (GetWidth() != 4 || GetHeight() != 4)
		throw "DoubleMatrix: Invalid DCT dimension";

	// X = At * Y * A;
	DoubleMatrix result = At;
	result *= *this;
	result *= A;

	*this = result;

	return *this;
}

std::ostream& operator<<(std::ostream& stream, const DoubleMatrix& m)
{
	unsigned int i, j;

	for (i = 0; i < m.GetHeight(); i++)
	{
		for(j = 0; j < m.GetWidth(); j++)
			stream << fixed << right << setw(8) << setprecision(3) << m.GetElement(i,j) << " ";

		stream << endl;
	}

	return stream;
}
